function[IRF,e,nwvar_IRF,sd_IRF]=LP_model(lhs,rhs,hor,HAC_sd_binary)
% Viet Hoang Dinh
% Monash University
% Edited December 2024

%%
N=size(rhs,2);
Phi=(rhs'*rhs)\(rhs'*lhs);
[msg, id] = lastwarn;
warning('off', id);
e=lhs-rhs*Phi;
if HAC_sd_binary==0
    nwvar=0;
elseif HAC_sd_binary==1
 nobs=size(lhs,1);
     HAC = 0;
       for l = 0:hor   %This portion of the loop is to calculate the autocorrelation for the Newey-West std error
          w_l = 1-l/(hor+1);
            for t = l+1:nobs
             if (l==0)   % This calculates the S_0 portion
                HAC = HAC  + e(t) ^2 * rhs(t,:)'*rhs(t,:);
             else        % This calculates the off-diagonal terms
                HAC = HAC + w_l * e(t) * e(t-l)* ...
                    (rhs(t,:)' * rhs(t-l,:) + rhs(t-l,:)' * rhs(t,:));
             end
            end
      end
            HAC = (1/(nobs-N)).*HAC;            
nwvar= diag(nobs.*((rhs'*rhs)\HAC/(rhs'*rhs)));

end
nwvar_IRF=nwvar(1); 
sd_IRF=sqrt(nwvar_IRF); %HAC-Newey-West standard errors for IRFs by Local Projections
IRF=Phi(1);             %IRFs by Local Projections
end